Kohn-Sham Density Functional Inspired Approach to Nuclear Binding 
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A non-relativisitic nuclear density functional theory is constructed, not as done most of the 
time, from an effective density dependent nucleon-nucleon force but directly introducing in the 
functional results from microscopic nuclear and neutron matter Bruckner G-matrix calculations at 
various densities. A purely phenomenological finite range part to account for surface properties 
is added. The striking result is that only four to five adjustable parameters, spin-orbit included, 
suffice to reproduce nuclear binding energies and radii with the same quality as obtained with the 
most performant effective forces. In this pilot work, for the pairing correlations, simply a density 
dependent zero range force is adopted from the literature. Possible future extensions of this approach 
are pointed out. 
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It is common use in nuclear physics to describe proper- 
ties of nuclei in mean-field approximation, using effective 
density dependent nucleon-nucleon interactions. Proto- 
types of these interactions are the so-called Skyrme-forces 
[l| which are of zero range and density dependent. There 
also exist finite range versions, like the Gogny force [2| 
and the relativistic mean field (RMF) approach [3|. The 
number of parameters which enter these effective forces is 
typically around ten and they are adjusted to reproduce 
finite nuclei and some nuclear matter properties. How- 
ever, recently also data from a theoretically determined 
neutron matter Equation of State (EOS) [4[ have been 
used as input (see also an older attempt in this direction 
in Q). Generally these forces give rise to an effective 
nucleon mass m* < m with typically m* /m ~ 0.7 in 
the non-relativistic framework. With these ingredients 
nuclear mean field theories are very successfull to de- 
scribe nuclear properties as, e.g. binding energies, radii, 
but also excited states, fission barriers, and many things 
more 0, 0] . These effective forces are still under debate 
and constantly refined to account for the ever increas- 
ing set of nuclear data 0, 0] ■ For instance the evolution 
with isospin, to be considered for nuclei approaching the 
neutron drip as involved in stellar nucleo-synthcsis, is an 
active line of research. Also, usually ( apart from the 
Gogny force), additional parameters are needed to char- 
acterise nuclear pairing. 

The main objective of the present study is to show that 
within an approach inspired by the Kohn-Sham Den- 
sity Functional Theory (KS-DFT), like it is presently 
still used extensively in, e.g^, condensed matter, chem- 
istry, and atomic physics [3, [l^, [ll| ( see, however, 
also [H, E 0] for updated versions of the universal 
exchange-correlation energy including gradients and ki- 
netic energy density dependences), a very efficient energy 



density functional for finite nuclei can be built up us- 
ing explicitly a fully microscopic input from nuclear and 
neutron matter calculations. We are, however, aware 
of the fact that the applicability of the KS-DFT ap- 
proach to self-bound systems, as nuclei, is not obvious 
and that presently vivid debate about this topic is go- 



ing on 



m EE Ea El, Ei- 



In this work we shall not 



contribute anything fundamental to this discussion but 
rather adopt a pragmatic point of view: try and see! 
The KS-DFT scheme is somewhat different from the 
usual Skyrme or Gogny approaches [H, Q where, as men- 
tioned, most of the time an effective force is constructed. 
However, several authors also interpret the Skyrme force 
as based on DFT [2(3, HH and a few authors have, in 
the past, constructed local Energy Density Functionals 
(EDF's) which cannot be related to underlying effective 
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for a review of differ- 
ent mean field calculations used in nuclear physics. In 
detail , however, all these nuclear approaches differ from 
what is practice in condensed matter physics (see below) . 
The basis of KS-DFT lies in the Hohenberg-Kohn (HK) 
theorem (26[, which states that for a Fermi system, with 
a non-degenerate ground state, the total energy can be 
expressed as a functional of the density p(r) only. Such 
a functional reaches its variational minimum when eval- 
uated with the exact ground state density. Furthermore, 
for practical reasons, in the KS-DFT method one intro- 
duces an auxiliary set of A orthonormal single particle 
wave functions ipi (r) , where A is the number of particles, 
and the density is assumed to be given by 
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p(r) = E ita , t \ipi(r,s,t)\ 



(1) 



where s and t stand for spin and iso-spin indices. The 
variational procedure to minimise the functional is per- 
formed in terms of the orbitals instead of the density. 
Usually in condensed matter and atomic physics the 
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HK functional E[p(r)] is split into two parts: E = 
T [p] + W[p] [8|. The first piece T corresponds to the 
uncorrelated part of the kinetic energy and within the 
KS method it is written as 

To = £E [d 3 r\VMr,s,t)\ 2 . (2) 

i,s : t J 

The other piece W[p] contains the potential energy as 
well as the correlated part of the kinetic energy. 

Then, upon variation, one gets a closed set of A 
Hartree-like equations with an effective potential, the 
functional derivative of W[p] with respect to the local 
density p(r). Since the latter depends on the density, 
and therefore on the ipi's, a self-consistent procedure is 
necessary. Still the equations are exact but they only can 
be of some use if a reliable approximation is found for the 
otherwise unknown density functional W[/>]. In the KS- 
DFT formalism the exact ground state wave function is 
actually not known, the density being the basic quantity. 

In nuclear physics, contrary to the situation in con- 
densed matter and atomic physics, the contribution of 
the spin-orbit interaction to the energy functional is very 
important. Non-local contributions have been included 
in DFT in several ways already long ago (see 0] for a 
recent review of this topic). Consequently, the spin-orbit 
part also can be split in an uncorrelated part E s o - plus a 
remainder. The form of the uncorrelated spin-orbit part 
is taken exactly as in the Skyrme [l| or Gogny forces • 

We thus write for the functional in the nuclear case 
E = T a + E s o - + E vnt + E c , where we explicitly split 
off the Coulomb energy Ec because it is a quite dis- 
tinct part in the Hamiltonian. It shall be treated, as 
usual, at lowest order, i.e. the direct term plus the ex- 
change contribution in the Slater approximation, that is 
Eg = (1/2) / / d 3 rd 3 r' p p (r)\r - t'\- 1 p p (t'), and Eg = 
-(3/4)(3/7r) 1 / 3 /d 3 r / 9 p (r) 4/3 with E c = Eg + Eg* and 
Pp/ n the proton/neutron density. 

Let us now discuss the nuclear energy functional contri- 
bution Ei nt [pn, Pp] which contains the nuclear potential 
energy as well as additional correlations. We shall split it 
in a finite range term E^[p n , p p ] to account for correct 
surface properties and a bulk correlation part E°^ t [p n , p p ] 
that we take from a microscopic infinite nuclear matter 
calculation e 29| as we will discuss below. Thus our final 
KS-DFT -like functional reads: 

E = T + E s -°- + E£ t + Ef n f + E C . (3) 

For the finite range term we make the simplest phe- 
nomenological ansatz possible 

Ef n f[p n ,P P ] = \Y,I [d 3 rd s r'p t (r)v ttt ,{r-r')p t ,(r') 

- I d3 m(r)pAr) (4) 

t,f J 



with t = proton/neutron and jt,t> the volume integral 
of ft, *'(?")■ The substraction in Jl]) is made in order not 
to contaminate the bulk part, determined from the mi- 
croscopic infinite matter calculation. Finite range terms 
have already been used earlier, generalizing usual Skyrme 
functionals (see e.g [HI EH HH ) • In this study, for the fi- 
nite range form factor Wt,t'( r ) we make a simple Gaussian 
ansatz: Kt.t'(r) = Vt.t'e~ r . We choose a minimum of 
three open parameters: V PiP = V n ,n = Vl,V UiP = V p . n = 
Vu, and r . 

The only undetermined and most important piece in 
© is then the bulk contribution E°^ t . In condensed mat- 
ter, chemistry and atomic physics, this quantity implies 
the EOS of interacting electrons and for the KS-DFT 
scheme, most of the time, it is obtained from Quantum 
Monte-Carlo (QMC) calculations, the results of which are 
then accurately represented by a fit function (not neces- 
sarily a polynomial, as we use later). As already men- 
tioned, we obtain E°% t from microscopic infinite matter 
calculations, using a realistic bare force, together with a 
converged hole- line expansion . We first reproduce by 
interpolating functions the correlation part of the ground 
state energy per particle of symmetric and pure neu- 
tron matters, and then make a quadratic interpolation 
for asymmetric matter. Finally the total correlation con- 
tribution to the energy functional in local density approx- 
imation reads: 

E%APp,Pn] = J d 3 r[P s (p)(l - /3 2 ) + P n (p)f3 2 ]p (5) 

where P s and P n are two interpolating polynomials for 
symmetric and pure neutron matter, respectively, at the 
density p — p p +p n , and (3 = (p n —p p )/ pis the asymmetry 
parameter. For P s we took ( x = p/po, with /?o=0.17 
fm -3 , see below ) 

Ps(p) = { 

{ P s (p ) + ax ■ (x - 1) + a 2 ■ (x ~ l) 2 x > 1 

(6) 

where the coefficients ( energy/particle in MeV) are given 
in Table I. The two forms match at x — 1 (p = p ) up to 
the second derivative. This functional form can be used 
up to p = 0.24 fm^ 3 , which is the interval where the 
independent fit of the microscopic calculation has been 
performed. 

A similar expression holds for P n , 

Pn(p)=I>i"V (7) 
k=l 

where again the coefficients are displayed in Table I, 
which is valid in the same density interval. The inter- 
polating polynomial for symmetric matter has been con- 
strained to allow a minimum exactly at the energy E j A 
= - 16 MeV and Fermi momentum kp = 1.36 fm _1 , i.e. 
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FIG. 1: EOS of symmetric and neutron matter obtained by 
the microscopic calculation (squares) and the corresponding 
polynomial fits (solid lines). 



Po— 0.17 fm 3 . This is within the uncertainty of the nu- 
merical microscopic calculations of the EOS 
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constrained fit was performed by keeping the EOS as 
smooth as possible, thus allowing for some very small 
deviations from the microscopic calculations below sat- 
uration density. An interpolating fit which goes exactly 
through the calculated EOS, as performed in [29], gives 
a not good enough saturation point (typically E/A = - 
15.6 MeV, kp=1.38 fm -1 ). In this context it should be 
noticed that the microscopic calculations of Ref. [32] have 
been modified, in order to get a good saturation point, by 
an ad hoc correction which amounts to more than 2 MeV 
in energy ( see,e.g. [33j for a presentation of the uncor- 
rected EOS of pH). As discussed in [33|], the low density 
behavior of the nuclear matter EOS is quite intricate and 
usually not reproduced by Skyrme and Gogny function- 
als ( see also ref. [HI), missing quite a substantial part 
of binding. We show our EOS for nuclear and neutron 
matter in Fig. 1. One sees that for all densities relevant 
to finite nuclei the EOS is very accurately reproduced by 
the interpolation (continuous lines). The bulk part E°% t 
of our functional, directly related to bare NN and NNN 
forces, is, therefore, determined once and for all and we 
will use it in (3) together with LDA. 

The only open parameters are, therefore, the ones con- 
tained in the finite range surface part, Eq.Q, and the 
strength of the spin-orbit contribution. We, thus, follow 
exactly the strategy employed in the above mentioned in- 
teracting electron systems. On the contrary, up to now, 
in nuclear physics, almost exclusively a different strategy 
was in use (see, however, Ref. [341] for a somewhat simi- 
lar approach to our): functionals like the one of ([3]), i.e. 
bulk, surface, etc., were globally parametrized with typi- 
cally of the order of ten parameters which, then, were de- 
termined fitting simultaneously some equilibrium nuclear 
matter (binding energy per particle, saturation density, 
incompressibiliy, etc.) and finite nuclei properties. How- 
ever, in this way, bulk and surface are not properly sepa- 



rated and early attempts used to miss important infinite 
nuclear matter properties, as, e.g. stability of neutron 
matter at high density [2j and other stability criteria. 
Modern Skyrme forces , like the Saclay-Lyon (SLy) ones, 
explicitly use the high density part {p/po > 0.65) of mi- 
croscopic neutron matter calculations for the EOS in the 
fitting procedure [1] and thus avoid collapse. Therefore, 
modern functionals usually reproduce reasonably well mi- 
croscopically determined EOS for neutron and nuclear 
matter (29| . Examples are, among others, the SLy- forces 
Q (see Fig.l in [2i|) and the Fayans functional (see e.g. 
Fig. 3 in @). It is thus evident that the procedure 
adopted in this work is different on a qualitative level 
from the usual and allows, via the fit, to reproduce very 
accurately the microscopic infinite matter results in the 
whole range of densities considered. This may be impor- 
tant for surface properties and neutron skins in exotic 
nuclei, what shall be investigated in the future. Let us 
also mention, in this context that for the electron systems 
mentioned before, the DFT practitioners demand to have 
extremly precise EOS as LDA input to their calculation. 

For open shell nuclei, we still have to add pairing. 
The formal generalization of the rigorous HK theorem to 
paired system has been given in Ref. [35| . In the present 
work our main objective is to set up the KS-DFT scheme 
for the non pairing part, thus we add pairing in a very 
simple way within the BCS approach. For this we simply 
take the density dependent delta force defined in Ref [36j 
for m = m* with the same parameters and in particu- 
lar with the same cutoff. As far as this amounts to a 
cutoff of ~ 10 MeV into the continuum for finite nuclei, 
we have to deal with single-particle energy levels lying in 
the continuum. We have simulated it by taking in the 
pairing window all the quasi-bound levels, i.e. the lev- 
els retainded by the centrifugal (neutrons) and centrifu- 
gal plus Coulomb (protons) barriers. This tratment of 
the continuum works properly, at least for nuclei not far 
from the stability valley as it has been extensively shown 
in [33] • A more sophysticated procedure, as, e.g, the 
one proposed by Bulgac et al [2J], may turn out to lead 
to better results. Since we do not want to concentrate 
on pairing here, this shall be investigated in forthcoming 
work. 

In our calculations the two-body center of mass cor- 
rection has been included in the self-consistent calcula- 
tion using the pocket formula, based on the harmonic os- 
cillator, derived in Ref438(] which nicely reproduces the 
exact correction as it has been shown in [32}. Our func- 
tional is now fully defined and, henceforth, we call it 
BCP-functional. 

In this exploratory investigation, we fitted two sets of 
parameters. We have considered only spherical nuclei 
which we chose as given below. The two fits were ob- 
tained in i) optimising ground state energies only, ii) op- 
timising ground state energies together with the charge 
rms radii. The open parameters, Vl, Vu, f*o of Eq.Q as 
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well as the spin-orbit strength Wo, are fitted to repro- 
duce the binding energies of the nuclei 16 O, 40 Ca, 48 Ca, 
56 Ni, 78 Ni, 90 Zr, 116 Sn, 124 Sn, 132 Sn, 208 Pb and 214 Pb 
in the case of the parameter set BCP1 and additionally 
the charge rms radii (r c = ^Jfp + 0.64 fm) of the nuclei 

16 0, 40 Ca, 48 Ca, 90 Zr, 116 Sn, 124 Sn, 208 Pb and 214 Pb for 
the parameter set BCP2. Experimental binding energies 
and charge radii are taken from References [40j and 41 1 
respectively. In table II we give the obtained parameter 
sets of fits (i) and (ii). 

Our intention is to compare our results from the BCP- 
functionals with the ones obtained from some of the 
most performant effective nucleon-nucleon forces avail- 
able, namely D1S 0], NL3 53 and SLy4 0. To this 
end, we have calculated the binding energies and charge 
radii of 161 even-even spherical nuclei (in line with the 
NL3 calculations reported in 53)- In Table III we re- 
port the energy and charge radii rms deviations between 
the corresponding experimental values and the theoret- 
ical ones obtained using the D1S force 5Hi the NL3 
parametrization 44 1, the Skyrme interaction SLy4 Q 
and our BCP1 and BCP2 functionals. We also display in 
Figs. 2-3 the differences between the theoretical and ex- 
perimental energies and charge radii in the range between 
16 Ne and 224 U calculated with the BCP1 functional in 
comparison with the same quantities obtained from the 
D1S effective force [iBj . 

In Table IV we display the two-neutron separation en- 
ergies for some magic nuclei predicted by the BCP1 and 
BCP2 functionals as compared with the same results pro- 
vided by D1S and NL3 as well as the corresponding ex- 
perimental values. 

We see that our theory is at least as performant as the 
other ones. In Table V we give the common bulk proper- 
ties of both BCP-functionals. Incompressibility modulus 
Koo and symmetry energy J plus its derivative L have 
quite acceptable values. For the average gap in 116 Sn we 
obtain 1.5 MeV. For the drip line nucleus 176 Sn we get 
E/A=-6.57 MeV compared to E/A=-6.U MeV (D1S) 
and E/A=-6.67 MeV (NL3). 

In a recent analysis, Bertsch et al.[46| have pointed out 
that also in Skyrme forces implicitly only four or five pa- 
rameters are relevant. For example, in the SLy forces [4J, 
six out of the twelve parameters are strongly constrained 
by nuclear and neutron matter properties. However, the 
philosophy in the present work is different and, as we said 
before, is analogous to what is done in condensed matter 
physics. That is we use as an input an analytic repre- 
sentation of a very accurate ab initio calculation of the 
nuclear/neutron EOS in the whole relevant range of den- 
sities, feed this into the functional via LDA and fit the 
remaining four parameters, which characterize the sur- 
face and spin-orbit contributions to the energy density, 
to finite nuclei properties. It is remarkable and, to some 
extend unexpected, that this scheme works indeed very 



well. It seems to give some credibility to the applicability 
of the KS-DFT scheme, even for self bound systems, as 
treated here (see also a discussion of this point in [181]). 



TABLE I: Coefficients of the polynomial fits P s , Eq. (0 and 
P n , Eq. to the EOS of symmetric and neutron matters. 



k 


bi s) (MeV) b^'(MeV) 


a fc (MeV) 


1 


-105.640069 


-43.985736 


-15.3563461 


2 


167.700968 


49.784439 


16.4197441 


3 


-181.762432 


-42.400650 


0.0 


4 


103.166047 


21.894382 


0.0 


5 


-22.4990207 


-4.3071179 


0.0 



TABLE II: Parameters of the Gaussian form factors and spin- 
orbit strength. 

r (fm) V L (MeV) Vu (MeV) W (MeV) 



BCP1 
BCP2 



1.05 
1.25 



-93.520 
-33.700 



-60.577 
-32.483 



113.829 
110.812 



TABLE III: Energies (in MeV) and charge radii (in fm) rms 
deviations. The numerical values of energies and charge radii 
calculated with the BCP1 and BCP2 functionals used to ob- 
tain the rms values of this Table are given in [471 ] 



BCP1 BCP2 PIS NL3 SLy4 
2.057 2.414 3.582 1.711 
0.028 0.020 0.020 0.024 



rmse 1.775 
rmsfl 0.031 



In conclusion, we have succesfully applied a Kohn- 
Sham Density Functional Theory inspired approach to 
ground state properties of nuclei, quite in analogy to 
what is usually done in e.g. condensed matter physics 
and for atoms and/or molecules [t|[l(|, however, at vari- 
ance with the procedure hitherto almost exclusively used 
in nuclear physics. This inspite of the fact that the ap- 
plicability of KS-DFT for self-bound systems is, so far, 
an unsettled problem. The bulk part of the functional is 
given for all relevant densities and asymmetries, once and 
for all, from microscopic results using the converged hole- 
line expansion based on realistic bare forces . To this 
a phenomenological surface part with three adjustable 
parameters is added. A fourth parameter is the strength 
of the spin-orbit interaction for which we get values close 
to the usual ones in Gogny and Skyrme forces. The mi- 
croscopic nuclear matter EOS has been fine tuned to pass 
through E/A = -16.00 MeV at saturation. This may be 
considered as the fifth adjusted parameter. In this pi- 
lot study we took for the pairing part of the functional 
the simplest possible procedure and adopted a previously 
adjusted density dependent d- force (3(|. The fit of open 
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TABLE IV: Two-neutron separation energies S-jn (in MeV) 
of some magic nuclei computed with BCP energy functionals, 
D1S force and NL3 parametrization Experimental values 
are also displayed. 





ICQ 


40 Ca 


48 Ca 


132 Sn 


208p b 


BCP1 


25.57 


26.62 


15.69 


11.94 


14.77 


BCP2 


25.82 


26.93 


15.97 


12.05 


14.79 


D1S 


29.57 


29.60 


15.68 


12.79 


14.89 


NL3 


28.41 


29.72 


16.36 


12.55 


13.99 


exp 


28.87 


28.93 


17.22 


12.56 


14.11 




I I I I I I I I I I I I I I I I I I I I I I I I I 

20 40 60 80 100 120 140 160 180 200 220 240 
A 



TABLE V: Common infinite nuclear matter properties of both 
BCP functionals 

B/A(MeV) po (fm~ 3 ) m/m* J (MeV) L(MeV) ^(MeV) 
^ilTOO 0T7 L00 33755 56739 2497 



parameters was done only for spherical nuclei. In view of 
the surprisingly small number of adjustable parameters, 
the results are very encouraging and well compete with 
the most performant mean field theories presently in use. 
In addition we have preliminary results which show that 
our BCP-functional also yields excellent results in the 
deformed case, again in very close agreement with those 
obtained with the Gogny D1S force |48j. Our study thus 
seems to support the applicability of the KS-DFT scheme 
also of self-bound systems. A more elaborate parameter 
search with different flexibilites than in Q is an impor- 
tant task for the near future. 

We are indebted to C.Maieron for valuable discus- 
sions. One of us (P.S.) is greatful to M.Casida for use- 
ful informations. This work is partially supported by 
IN2P3-CICYT and INFN-CICYT. X.V. acknowledges 
Grant Nos.FIS2005-03142 (MEC, Spain and FEDER) 
and 2005SGR-00343 (Generalitat de Catalunya). 



FIG. 2: Differences between the theoretical and experimental 
energies. Calculations are performed using the BCP1 func- 
tional (open circles) and the D1S parameter set (crosses) [45| . 
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FIG. 3: Differences between the theoretical and experimen- 
tal charge rms radii. Calculations are performed using the 
BCP1 energy functional (open circles) and D1S parameter 
set (crosses) [iq] . 
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